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CHAPTER  I:  INTRODUCTION 


Purpose 

The  purpose  of  this  research  is  to  establish  a  preliminary  model  to 
investigate  the  crack  propagation  history  in  a  given  pavement  system.  This 
is  accomplished  through  a  multi-step  process.  First,  a  suitable  program  is 
used  to  calculate  the  stress  distribution  in  the  pavement  layers  under  the 
load,  in  this  study  an  elastic  layer  program.  Second,  through  successive 
application  of  a  finite  element  program,  the  stress  intensity  factor  as  a 
function  of  crack  length  is  determined.  Third,  using  the  stress  intensity 
factor  distribution,  the  number  of  load  cycles  required  to  advance  the  crack 
a  given  increment  is  calculated.  These  increments  carry  the  crack  from  its 
initial  to  its  final  value,  which  defines  failure.  The  load  cycles  are 
calculated  using  the  Paris  equation.  The  remainder  of  the  report  is  taken  up 
by  a  detailed  description  of  the  problem  and  an  in-depth  account  of  the 
solution  process  and  results. 

Pavement  System 


The  pavement  system  used  in  this  study  consists  of  a  combination  of  three 
layers:  a  subgrade,  a  base,  and  a  surface.  Each  layer  has  variable  values 
for  the  modulus  of  elasticity  and  the  thickness.  This  information  is 
summarized  in  Table  1. 


The  subgrade  thickness  is  listed  as  50  inches.  This  is  an  arbitrary 
value  selected  so  that  the  subgrade  thickness  does  not  affect  the  solution. 
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Table 

1.  Pavenent  Data  Used 

in  Developing  Structural  Data. 

LAYER 

MODULUS  VALUES. PSI 

THICKNESSES, IN. 

POISSON-S  RATIO 

ovcrl ay 

3S0000.0 

0 

3 

6 

0.3S 

Ifc, 

base 

400000.0 

8 

0.3S 

800000.0 

12 

1200000.0 

16 

subgrade 

4000.0 

SO 

0.3S 

9000.0 

MODULUS 

CONDITIONS 

MODULUS 

CONDITION  BASE  MODULUS. PSI  SUBGRADE 

MODULUS. PSI 

: 

1 

1200000.0 

4000.0 

2 

1200000.0 

9000.0 

3 

800000.0 

4000.0 

4 

800000.0 

9000.0 

5 

400000.0 

4000.0 

6 

400000.0 

9000.0 

m 


In  the  analysis  of  the  pavement  system,  all  possible  combinations  of  the 
values  in  Table  1  are  used  to  develop  the  influence  of  pavement  properties  on 
fracture , 


Loading 

The  loading  shown  in  Figure  1  approximates  that  of  an  F-4  aircraft.  It 
consists  of  a  total  load  of  27,000  pounds.  This  load  is  distributed  over  a 
circle  with  a  radius  of  5.7  inches  to  provide  a  pressure  of  265  psi. 


Solution  Philosophy 


Strictly  speaking,  this  problem  is  a  three  dimensional  axi-symmetr ic 
one.  However,  because  crack  propagation  itself  is  generally  a  plane  strain 
phenomenon,  it  is  considered  preferable  to  rely  on  a  simpler  solution 
scheme.  The  approach  here  is  to  use  a  two-dimensional  plane  strain  finite 
element  program  to  obtain  approximate  stress  intensity  factors  developing 
under  the  stresses  calculated  from  an  elastic  layer  analysis.  for  a  single 
load,  as  assumed  here,  the  results  are  satisfactory. 

The  cornerstone  of  the  solution  procedure  is  the  computation  of  the 
stress  intensity  factor  as  a  function  of  crack  length.  This  is  done  with  the 
finite  element  program  described  in  a  later  chapter  and  with  the  aid  of 
superposition  techniques.  When  dealing  with  boundary  values  over  an  infinite 
region  such  as  a  pavement  structure,  the  problem  can  be  analyzed  as  the 
superposition  of  two  problems  (^)  as  illustrated  in  Figure  2.  In  the  first 
phase,  one  applies  the  total  loading  to  the  pavement  with  no  crack.  For  the 
second  phase,  the  crack  is  included  in  the  pavement.  The  stresses  where  the 


Surface 


I 


Base 
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P=27,000  pounds 
a=265  psi 
r>=5.7  Inches 


Figure  1.  Loading  Condition  Modelled  With  Finite  Element 
Fracture  Program. 


crack  would  be  in  the  first  phase  are  applied  in  the  reverse  direction  to  the 
crack  in  the  second  phase.  The  stress  intensity  factor  found  from  the  second 
phase,  although  of  opposite  sign,  is  equal  in  magnitude  to  that  of  the 
overall  problem.  Thus  the  task  of  finding  Kj  becomes  one  of  discovering 
the  correct  stress  distribution  due  to  loading  at  the  crack,  applying  the 
reverse  stresses  to  the  crack,  and  analyzing  the  pavement. 

A  further  difficulty  arises  when  using  the  above  superposition.  As  will 
be  noted,  stresses  cannot  be  applied  to  the  portion  of  the  crack  which  lies 
within  the  crack  tip  finite  element.  Thus,  a  correction  factor  is  required 
to  account  for  this.  The  correction  factor  is  given  by  (1)  as: 


C 


k 


IT 


a  (§)  di 
e 

/T 


(1) 


0 

where  §  is  the  distance  away  from  the  crack  tip  and  o^(§)  is  the  stress  to  be 

applied  to  the  portion  of  the  crack  that  lies  within  the  crack  element.  If  z 

is  sufficiently  small,  a  (§)  can  be  assumed  to  be  a  constant,  and  equation 

e 

(1)  reduces  to  (1): 


(2) 


The  final,  correct  stress  intensity  factor  is  found  from  Cjj  and  Kj 
computed  in  superposition  phase  two  by 


final  ”  computed 


(3) 


Crack  propagation  is  calculated  using  this  stress  intensity  factor  in 


numerical  integration  of  the  Paris  equation.  The  Paris  equation  consists  of 
the  following: 


dc/dN  -  AKj" 


(A) 


where 

dc  =  differential  crack  extension 

dN  *  differential  increase  in  the  number  of  load  cycles 
Kj  *  stress  intensity  factor 

A,n  =  dimensionless  material  constants  determined  experimentally 

Manipulation  of  the  Paris  equation  gives: 


dN  •  dc/AKi" 


(5) 


from  which  ANf,  the  number  of  load  cycles  required  to  advance  the  crack  an 
increment  can  be  calculated  as: 


o 

Writing  the  integral  in  terms  of  c/b  leaves: 


AN 


f 


d(c/b) 
A  K," 


(6) 


(7) 


When  the  current  crack  length  equals  the  crack  length  at  failure,  the  final 
value  of  Nf  is  reached. 
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The  total  number  of  load  cycles  required  to  advance  the  crack  from  the 
initial  crack  length  to  the  current  crack  length,  Nfj,,  is  found  by  summing 
the  ANf  values  for  the  n  crack  increments  which  total  the  current  crack 
length.  Thus, 


=  E  an 


,  *  N,.  +  AN, 

fn  f(n-l)  fn 


Note  that  Kj  is  considered  a  function  of  (c/b)  now,  not  crack  length.  This 
allows  regression  equations  to  be  developed  for  the  range  of  materials  and 
loading  conditions  which  lessens  the  reliance  on  the  finite  element  computer 
code . 

In  this  analysis  the  crack  is  assumed  to  form  directly  beneath  the  center 
line  of  the  load,  on  the  line  of  symmetry.  It  is  assumed  to  originate  at  the 
interface  of  the  base  and  the  subgrade  and  to  propagate  upward  through  the 
base.  The  greatest  crack  length  at  which  the  pavement  is  assumed  to  have 
failed  is  the  full  thickness  of  the  base.  If  c  ■  crack  length  in  the  base 
and  b  =  thickness  of  the  base,  then  c/b  «  1.0  is  the  maximum  value  at  which 
failure  can  occur.  This  would  lead  one  to  believe  that  the  life  of  the 
pavement  ranges  between  c/b  values  of  0.0  and  1.0.  However,  this  is  true 
only  if  two  conditions  exist.  First,  Kj,  the  stress  intensity  factor,  must 
be  positive  for  all  c/b.  Second,  the  critical  value  of  the  stress  intensity 
factor,  Kj(;,  must  not  be  exceeded.  In  the  first  case,  if  Kj  is  negative, 
then  the  crack  cannot  propagate  in  that  region  and  that  region  cannot  figure 
in  the  fracture  life  of  the  pavement.  The  areas  defining  the  fracture  life 
of  the  pavement  when  K  is  negative  are  illustrated  in  Figure  3,  where  typical 
Kj  distributions  for  this  problem  are  depicted.  Case  I  and  II  are 
self-explanatory.  In  case  III,  it  is  assumed  that  some  perturbation  advances 


«•.  r.  / 


the  crack  to  (c/b)^,  so  that  the  crack  can  propagate.  Here  (c/b)^  is 
considered  small  so  that  this  is  readily  possible  and  handled  with  the 
concept  of  starter  flaws.  If  (c/b)2  is  not  small,  then  the  crack  may  never 
propagate.  In  case  IV,  the  crack  propagates  to  (c/b)^  in  the  normal 

fashion  and  once  there  encounters  negative  Kj  values  and  stops.  It  is 
unlikely  that  the  crack  will  somehow  be  advanced  suddenly  through  the 

relatively  extensive  negative  Kj  region,  and  if  it  is,  the  pavement  will  be 
nearly  destroyed  and  can  be  said  to  have  failed.  Thus  (c/b)^  is  the  upper 
limit  of  the  fatigue  life  in  this  case.  One  should  note  that  if  (c/b)^  is 
significantly  less  than  1.0,  then  exhausting  the  fatigue  life  may  not 
constitute  failure  of  the  pavement.  The  pavement  may  still  maintain  its 
structural  integrity. 

If  Kj(’  is  exceeded,  then  the  crack  instantly  propagates  to  a  point 

where  Kj  ■  Kj(;.  The  fatigue  life  would  then  range  between  c/b  values 

where  Kj  is  less  than  or  equal  to  Kj^.  If  Kj  never  falls  below  Kjg 

again,  then  the  crack  instantly  propagates  to  failure.  Unfortunately,  Kjc 
values  for  the  pavements  in  this  problem  are  unknown,  so  there  is  no  way  to 
incorporate  them  into  the  solution  algorithm.  Therefore,  for  purposes  of 
analysis,  it  is  assumed  that  is  greater  than  Kj  for  all  c/b. 


Solution  Sequence 


the  solution  sequence  briefly  described  here  will  be  described  in  detail 
in  the  following  chapters.  The  stress  distribution  calculations  will  not  be 
described,  as  this  program  can  be  replaced  with  any  program  which  calculates 
stresses  and  strains  under  a  loading.  Other  elastic  theory  programs  can  be 
used  which  will  allow  multiple  wheel  loadings  to  be  used  to  produce  the 


stress  distribution  which  is  then  used  in  the  finite  element  fracture 
program.  The  Hybrid  crack  tip  element  of  Pian  and  Tong  used  in  a  plane 
strain  finite  element  program  is  used  to  determine  the  stress  intensity 
factors  resulting  when  a  cracked  base  course  is  subjected  to  the  calculated 
stress  distribution.  Regression  equations  have  been  developed  in  this  study 
to  reduce  the  use  of  the  finite  element  program  for  simple  pavement 
structures.  Finally,  the  Paris  equation  is  used  in  a  program  to  calculate 
the  number  of  loads  to  failure  under  the  given  stress  intensity  distribution. 
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CHAPTER  II:  THE  FINITE  ELEMENT  PROGRAM 


The  Singularity  Problem 


For  linearly  elastic  plane  strain  and  plane  stress  problems,  in  the 
vicinity  of  a  crack  tip  the  stress  varies  as  1/  /r,  where  r  is  the  radial 
coordinate  of  any  point  in  the  plane  measured  from  the  crack  tip.  Thus,  at 
the  crack  tip  where  r  «  0  the  stress  is  singular.  In  the  equations  for 
stress,  the  coefficient  of  the  singular  1/  /T"  term  is  called  the  stress 
intensity  factor  and  is  representative  of  the  "strength"  of  the  singularity. 
It  has  units  of  (force/area)*  /length.  There  are  two  main  types  of  crack 
which  govern  the  behavior  of  linearly  elastic  plane  stress/strain  problems. 
These  are: 

1)  mode  I  or  opening  mode 

2)  mode  II  or  in-plane  shearing 
These  modes  are  illustrated  in  Figure  4. 

For  isotropic  materials,  the  singular  terms  of  the  stress  distribution  in 
the  vicinity  of  a  crack  tip  are  (_3) : 
mode  I: 
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In  many  formulations,  the  coefficient  involving  K  in  front  of  the  vectors 
containing  trigonometric  terms  in  equations  (8)  and  (9)  is  given  as  (K//2r). 
The  finite  element  program  used  in  this  research  uses  equations  (8)  and  (9) 
as  is.  With  the  use  of  these  equations  it  is  necessary  to  multiply  the 
results  of  the  finite  element  program  by  /tT  in  order  to  match  any  results 
obtained  from  stress  distributions  based  on  the  alternative  coefficient.  It 

is  a  minor  point  but  mix-up  in  definitions  can  lead  to  confusion  and 

incomprehensible  results. 

The  stress  intensity  factors  Kj  and  Kjj  are  directly  proportional  to 

the  magnitude  of  applied  loading  and  are  dependent  on  the  geometry  of  the 

structure,  the  size  and  shape  of  the  crack,  and  the  nature  of  the  applied 
load ing . 

Displacement  fields  in  the  vicinity  of  a  crack  tip  are  also  functions  of 
the  stress  intensity  factors  and  can  be  written  as  : 
mode  I: 


>/8G 


mode  II: 


r  (2K  -  l)cose/2)  -  cos(3e/2)  -I 


(2K  ♦  l)sin(e/2)  -  sin(3e/2)-' 


r  (2K  +  3)sin(e/2)  +  sin(3e/2)  -t 


(10) 


<>41  ^  )/8G 


(11) 


-(2k  -  3)cos(6/2)  -  cos(30/2)*' 

In  equations  (10)  and  (11),  u  and  v  are  the  displacements  along  the  x  and 
y  axes  in  Figure  I-l ,  G  is  the  shear  modulus,  k  ■  (3  -  4v )  for  plane  strain 


and  (3  +  '^)/(l  +  v)  for  plane  stress,  and  v  is  Poisson's  ratio. 


The  stress  intensity  factor  is  also  related  to  the  strain  energy  release 
rate,  that  is  the  change  in  strain  energy  in  the  structure  per  unit  distance 
of  crack  extension.  For  linearly  elastic  plane  strain/stress  problems  the 
strain  energy  release  rate  is  given  by  O) : 
mode  I : 

£  I  «  +  1 )/8G  (12) 

mode  II: 

^II  *  K^ijCk  1  )/8G  (13) 

where  k  and  G  are  as  defined  previously. 

The  computer  program  provides  for  the  calculation  of  the  stress  intensity 
factor  in  two  ways.  The  first  and  far  more  emphasized  way  is  to  utilize  eqs. 
(8)-(ll)  and  solve  for  the  stress  intensity  factor  directly.  The  second 
method  is  to  run  the  program  twice,  the  second  time  with  a  slightly  longer 
crack,  note  the  difference  in  strain  energy,  and  apply  equation  (12)  and 
(13).  This  option  is  there  if  the  analyst  wants  to  take  advantage  of  it,  but 
it  is  not  the  main  thrust  of  the  program  and  is  limited  in  that  if  both  mode 
I  and  mode  II  are  present  equations  (12)  and  (13)  are  not  applicable. 

The  justification  for  using  elements  with  assumed  stress  and  displacement 
distributions  that  have  the  1/  /r*  singularity  built  in  is  based  on  the 
convergence  rate  of  problems  with  singularities.  The  convergence  of  such 
problems  has  been  shown  to  be  of  order  h,  where  h  is  the  maximum  size  of  the 
elements  used  in  the  solution.  The  convergence  rate  is  independent  of  p,  the 
order  of  the  complete  polynomial  used  in  the  interpolation  functions  for 
stress  and  displacement.  The  quantities  for  which  the  convergence  rate  was 
established  are  the  strain  energy  U  and  the  stress  intensity  factors  Kj  and 
Kii.  Given  the  above,  in  order  to  achieve  good  results  for  U  or  K  using 
elements  without  singularities  built  in,  the  number  of  elements  needed 
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becomes  impractical  because  results  for  relatively  large  h  are  highly 
inaccurate.  This  is  the  value  of  solutions  using  elements  with  singularities 
in  the  region  around  the  crack  tip.  Although  they  converge  no  faster,  for 
relatively  large  element  sizes  they  give  acceptable  results.  Indeed,  making 
the  elements  that  include  singularities  too  small  can  have  an  adverse  effect 
on  the  results.  This  peculiarity  occurs  for  this  reason.  If  the  elements 
containing  singularities  are  small  enough,  the  region  in  the  structure 
significantly  affected  by  the  crack  will  extend  beyond  the  scope  of  the 
singular  elements  and  cause  errors  in  adjacent  non-singular  elements.  In  the 
limit  of  mesh  refinement,  as  all  element  sizes  go  to  zero,  the  beneficial 
effects  of  elements  with  singularities  disappear,  leaving  only  a  solution 
based  on  conventional  elements. 


The  Crack  Element 


The  crack  element  first  developed  by  Pian  and  Tong  (^) ,  and  which  is 
shown  in  Figure  5  is  the  element  used  in  the  finite  element  program.  It  must 
be  of  sufficient  size  to  insure  good  results,  as  mentioned  above.  The  crack 
element  is  developed  using  a  functional  and  methods  which  are  a  combination 
of  both  assumed  stress  and  displacement  models.  As  such,  in  this 
development,  internal  stresses  and  displacements  as  well  as  boundary 
displacements  and  tractions  can  be  assumed  independently.  The  necessary 
equilibrium  and  compatibility  requirements  are  incorporated  into  the  Euler 
equations,  which  result  from  the  stationary  condition  of  the  functional.  The 
Euler  equations  can  be  satisfied  either  exactly  or  in  a  average  sense,  the 
latter  condition  being  the  source  of  approximation  in  the  finite  element 
method . 
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In  the  previous  section  a  term  is  included  in  the  equations  for 
displacements.  It  is  not  desirable  to  introduce  this  term  into  the  assumed 
boundary  displacements,  for  two  reasons.  First,  it  can  be  incorporated  into 
the  assumed  internal  displacements.  Second,  one  of  the  virtues  of  this  crack 
element  is  that  it  can  be  used  in  conjunction  with  conventional  finite 
elements,  and  if  a  /T  term  were  included  in  its  boundary  displacement, 
boundary  compatibility  between  it  and  a  conventional  element  could  not  be 
insured.  Of  course,  insuring  this  compatibility  causes  internal  and  boundary 
displacements  for  the  crack  element  to  be  incompatible,  a  source  of  error. 

The  functional  used  in  this  development  is,  for  plane  problems. 


TT 

m 


f 


(  U^-U^)T^ds- 


3A 


m 


where 


KH 


A 

m 


]dA  (14) 


A„  =  area  of  the  m^h  element 
3  A„  *  the  boundary  of 

Sn,  *  the  portion  ofSAj,  over  %diich  tractions  are  prescribed 

Sy  =  the  portion  of3An,  over  which  displacements  are  prescribed 

Cijki  =  the  elastic  compliance  tensor 

u£  =  a  smooth  internal  displacement  field 

U£  =  assumed  boundary  displacements 

u£  *  prescribed  boundary  displacements 

T£  ■  assumed  boundary  tractions 

T£  *  prescribed  boundary  tractions 

U£  is  assumed  such  that  u£  equals  u£  on  •  Note  that  the 

thickness  t  is  absent  from  equation  (14)  as  unity  is  assumed,  but  can  be 
multiplied  in  if  the  thickness  is  variable  variable. 


(15a) 


Euler  equations  for  the  functional 

in  equation  (7)  are 

in  Am 

in  Am 

on  Am 

Ui=Ui 

on  Am 

T.=T. 

on  m 

(15b) 

(15c) 

(15d) 

(15e) 


1  1 

Note  that  the  body  force  is  excluded  for  simplicity,  and  that  vj  is  a 
direction  cosine. 

Equation  (15a)  represents  compliance  with  stress-strain  laws  by  assumed 
stresses  and  displacements,  equation  (15b)  represents  the  satisfaction  of 
internal  equilibrium,  equation  (15c)  represents  compatibility  of  assumed 
stresses  and  boundary  tractions,  equation  (15d)  represents  compatibility  of 
internal  and  boundary  displacements,  and  equation  (15e)  represents 
satisfaction  of  boundary  conditions  by  boundary  tractions.  The  crack  element 
provides  for  internal  equilibrium  and  compatibility,  interelement 
compatibility,  and  boundary  conditions  on  displacements  identically.  The 
exact  solution  satisfies  all  the  Euler  equations.  Also,  if  the  solution  is 
exact,  then  the  stationary  condition  on  equation  (14)  will  provide  for 
interelement  equilibrium  and  compatibility.  Otherwise,  interelement 
equilibrium  is  satisfied  only  in  a  work  equivalent  sense,  while  interelement 
compatibility  will  be  satisfied  by  assuming  the  same  boundary  displacement 
functions  for  all  elements,  crack  or  conventional. 

The  element  stiffness  matrix,  whose  development  is  discussed  in  the 


following  pages,  is  derived  using  complex  variable  techniques. 
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To  begin  with,  let  z=x  +  iy,  where  x  and  y  are  as  depicted  in  Figure  4. 
As  shown  previously,  (a)~  //F  &  (^)  ,  so  that  stresses  and  displacements 

can  be  expressed  as  O ) : 

ay+o^=2[(i)' (2)+<t>' (z)  ]  /  (16) 

Oy-a^+ Z.  0x^=2  [Z<|>"(Z)+x'(Z)]  t 

Z^(U+iV)=ri4)(Z)-Z$'  (Z)-IhZ)  (17) 


where  =  E/2(l  +  v),  n =  3-4  v  for  plane  strain  and  (3  -  v)/(l  +  v)  for  plane 
stress.  E  and  v  are  Young's  modulus  and  Poisson's  ratio,  respectively.  It 


can  be  seen  that  y  ^  G  and 


T)  =  k  in  Section  I.  (  )'  denotes 


differentiation  and  (  )  denotes  the  complex  conjugate.  (f)  and  ij;  are  analytic 
functions.  The  above  definitions  of  stress  and  displacement  satisfy  Euler 
equations  (15a)  and  (15b).  In  addition,  boundary  tractions  are  chosen  in 
this  development  so  that  they  satisfy  equation  (15c). 

In  order  to  choose  proper  stresses  and  displacements  for  the  crack 
element  that  account  for  singularities  of  all  order,  the  following  mapping 
function  is  introduced  (^) . 

Z  =  W  (£)  =  (18) 


£*  (19) 

with  -Tr/2  <  arg  /2  and  -  X  argz  < On  the  £  plane,  the  element  lies 

on  the  region  where  the  real  part  of  £  is  positive  and  the  crack  lies  on  the 


imaginary  axis. 


and  }p  are  analytic  functions  of  £,  enabling  simple 


polynomials  in  terms  of  £  to  be  used  in  the  finite  element  solution. 


Using  the  mapping  given  above,  equations  (16)  and  (17)  become 


a  +a  -4R  [0'C£)/w'(£)l 

X  y  e 

a  +0-2  a  -2  (£)/w' (£)]’+l^•  (£)}/w’ (£)  (20) 

X  y  i 

2u(+iv)=n<|)(X)-wCt)a'  (£)/w' 

In  order  to  satisfy  the  stress  free  condition  on  the  crack  tip  given  by 

(3). 

a<£)-hv(£)VW)7^P^  +  T^^0  (21) 

the  following  form  of  is  chosen 

<P  '  (£) /w’  (£)  (22 ) 


By  using  equations  (15c),  (19),  (20),  and  (22)  all  of  the  Euler  equations 
except  for  equation  (15d)  are  satisfied  by  this  crack  element  model. 
Substituting  equations  (21)  into  (14)  gives  in  matrix  form  (_4) : 


TT  «  /  (T)'’^(U)ds- 


in  which 


/  (T)’^(U)d8 


r  0  V  +0  V  T 

1  ^  1 

XX  xy  y  1 

Lt  J 

L  0  V  +0  V 

2 

xy  X  y  y 

boundary  tractions 


=  internal  displacements 


boundary  displacements 


For  the  derivation  of  the  element  stiffness  matrix,  the  following  forms 
of  0(L)  and  ’^P  (£)  may  be  assumed;  (4) 


(25) 


N 

♦«>-  E  tr 

N  _  ..  si 

<)'(£)=- E  [b;(-i)^+i  bl]  ■ 
j=l  2 

where  N  is  a  finite  integer  and  bj  *  6j  +  i  Bjj+j  with  the  B's  being 
real  constants  and  ^n+2  “  This  is  because  <t>*  i  ^  and  0  give  no 

contribution  to  the  stresses  in  equations  (20). 

Using  equations  (20),  (24),  and  (25)  one  can  express  the  following;  (^) 

(  T  )  =  [  R  1(  B  )  (26) 

(  u  )  =  [  U  ](  B  ) 

where  (  B  )  includes  components  Bj,  B2, _ _  B2N  excluding  Bn+2" 

boundary  displacements  (  u  )  are  expressed  as 

(  u  )  »  [  L  ](  q)  (27) 

(  q  )  is  the  vector  of  nodal  displacements  and  [  L  ]  is  an  interpolation 
matrix  defined  on  3A^.  [  L  ]  is  such  that  boundary  displacements  for  the 

crack  element  and  adjacent  conventional  elements  are  the  same.  For  instance, 
if,  as  in  the  finite  element  program,  linear  boundary  displacements  are 
employed,  then  between  any  two  nodes  on  the  crack  element  [  L  ]  will  be  such 
that  (  u  )  has  a  linear  variation. 

The  number  of  nodes  in  the  crack  element  used  in  the  finite  element 
program  varies  from  5  to  9  depending  on  symmetry,  as  shown  in  Figure  5. 
Boundary  displacements  are  not  assumed  along  the  crack  edge  of  the  element 
because  tractions  are  zero  and  there  is  no  adjacent  element.  Because  of 
this,  in  the  program,  there  is  no  way  to  load  the  portion  of  the  crack  within 
the  element  and  there  is  no  way  to  control  the  crack's  displacement  within 
the  element. 

Substitution  of  equations  (26)  and  (27)  into  equation  (23)  gives 
’'m  -  (  B  )T[  G  ](  q  )  -  1/2  (  B  )T[  H  ](  B  )  (28) 


in  which 


[  G  ]  -/  [  R  )T[  L  ]ds  (29) 

■'sAm 

[  H  ]  -  1/2  /  ([  R  ]T[  u  ]  +  [  U  iTl  R  ])ds  (30) 

•'BAm 

As  the  B  's  can  be  assumed  independently  from  surrounding  elements,  the 
stationary  condition  of  equation  (28)  with  respect  to  (  B  )  gives 

[  H  ](  6)  -  [  G  ](  q  )  (31) 


or 

(  6)  “  [  H  J-M  G  ](  q  )  (32) 

With  equation  (31),  one  can  substitute  back  into  equation  (28)  and 
eliminate  (  B  ).  One  does  not  have  to.  Instead,  one  could  solve  for  (  B  )'s 
simultaneously  with  the  (  q  )'s.  However,  in  the  finite  element  program  the 
(  B  )'s  are  eliminated. 

Substituting  equation  (32)  into  equation  (28)  gives  (3^) 

=  l/2(  q  )T[  k  ](  q  )  (33) 

where  the  element  stiffness  matrix  for  the  crack  element  is 

(  k  ]  -  [  G  ]T[  H  G  ]  (34) 

After  global  assembly  and  solution  for  the  global  displacement  vector, 
the  stress  intensity  factors  can  be  found  from  the  now  known  (  q  ) .  The 
stress  intensity  factors  can  be  shown  to  be  related  to  (  B  )  in  the  following 
manner:  (4) 

Ki  =  (6i)/r 

(35) 

•^II  “  (  6n+1) 

Because  (  B  )  is  related  to  (  q  )  by  equation  (32),  one  can  write 

Ki  *  (  Bi)T(  q  ) 


*^II  “  (  ®II  9  ) 


m  ^  m  if.,,  A.i  1^4  mJK  -  p  ^  m.  r  *.  •:  «  < 


(36) 
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It  appears  Chat  (Bj)  is  the  first  row  of  [  H  ]  [  G  ]  multiplied 

by  and  that  (  is  the  (N+1)^^  row  of  [  H  ]“^[  G  ]  multiplied 

by  /2~ as  is  shown  in  the  program. 

The  reference  followed  in  Che  development  of  this  chapter  gives  details 
as  to  how  the  integrations  implied  by  equations  (29)  and  (30)  are  to  be 
accomplished  in  terms  of  complex  variables,  and  it  gives  certain  properties 
of  [  H  ]  (namely  that  its  upper  right  and  lower  left  quadrants  are  blocks  of 
zeros)  that  make  it  more  efficient  to  invert  (4), 

The  finite  element  program  uses  the  crack  element  of  the  previous  section 
in  conjunction  with  GST  triangular  elements  and  4-CST  quadrilateral  elements 
obtained  by  condensing  the  middle  node.  In  fact,  the  program  is  essentially 
Che  program  in  Desai  and  Abie's  finite  element  text  (  ^  )  modified  to 
incorporate  Che  crack  element. 

A  complete  discussion  of  the  entire  program  and  its  input  requirements  is 
given  in  Appendix  B. 

Operational  Parameters 

It  is  convenient  to  define  certain  parameters  which  are  used  to  describe 
and  discuss  the  crack  element.  The  first  parameter  is  c  . 

c  “(Length  of  crack  within  crack  element)/(i-otal  Crack  Length)* 

The  second  parameter  is  Che  ratio  a/1. 

a/1  “(Length  of  crack  within  crack  element)/(Length  of  Crack  Element) 

The  last  parameter  is  the  ratio  c/b. 


c/b  “(Total  crack  length) /(pepth  of  Base> 


The  element  shovm  in  Figure  6  is  used  for  all  the  finite  element  meshes 
in  the  solution  procedure  yet  to  be  described.  Because  of  the  symmetry  of 
the  problem,  the  five  node  element  is  used  here.  It  is  desirable  to  use  one 
element  for  all  the  meshes  so  that  one  can  create  and  build  meshes  quickly 
and  easily.  This  particular  element  is  chosen  for  various  reasons.  First, 
it  gives  a  reasonable  value  of  a/1  for  all  c/b  ratios.  A  plot  of  a/1  versus 
percent  error  in  Kj  is  shotm  in  Figure  7  for  the  crack  element  applied  to 
the  Bowie  crack  problem  •  The  above  element  has  a  constant  a/1  ratio  of 
0.5.  It  is  seen  that  this  a/1  ratio  yields  about  a  3  percent  error,  which  is 
relatively  small. 

A  second  reason  that  the  element  of  Figure  6  is  used  is  that  that  as  the 
crack  propagates,  e  remains  within  reasonable  limits  for  accuracy.  This  can 
be  seen  from  Figure  8,  %»hich  is  a  plot  of  e  versus  Kj  for  the  Bowie  crack 
problem  (_5) .  The  element  used  in  this  research  is  the  2  element.  From 
this  plot,  a  value  of  e  >  0.2  is  needed  to  insure  reasonable  accuracy  in 
Kj.  In  the  solution  process  for  the  problem  of  this  report,  a  maximum 
value  of  c/b  of  1/2  is  used,  and  a  minimum  c/b  value  of  1/12  is  used.  Also, 
the  greatest  b  value  is  16  inches,  and  the  least  is  8  inches.  Thus,  c,„ax 
is  8  inches,  and  c^jin  2/3  inches.  t  is  therefore  such  that 
0.03125  <  0.375,  and  within  the  range  of  acceptable  accuracy  shown  in 
Figure  8. 

The  above  two  reasons  for  using  the  element  shorn  in  Figure  6  are 
important,  but  many  elements  which  satisfy  the  a/1  and  criteria  could  be 
created.  The  main  reason  to  use  the  Figure  6  element  is  that  it  is  small 
enough  so  that  reverse  stresses  can  be  applied  to  enough  of  the  crack  length 
to  make  the  solution  process  viable.  For  instance, 


Cm,--  is  2/3  inches. 
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Figure  8.  Illustration  of  Accuracy  Related 
to  Crack  Tip  Length  (4) 


with  a  *  1/4  inches,  as  in  Figure  6,  one  can  apply  tractions  to  5/12  inches, 
or  62.5%,  of  the  crack  length.  Remember,  one  cannot  apply  tractions  to  the 
portion  of  the  crack  within  the  crack  element.  The  62.5%  value  can  be 
considered  to  be  a  significant  amount  and  meaningful  results  can  be  obtained 
from  it.  However,  as  a  increases,  this  percent  value  decreases  and  quickly 
becomes  insignificant.  If  a  =  c ,  or  1,  then  there  is  no  crack  length 
available  to  apply  reverse  stresses  to,  and  the  solution  philosophy  outlined 
earlier  1  is  impossible. 

There  will  be  five  values  of  c/b  used  in  the  finite  element  analysis  of 
each  combination  of  base,  overlay,  and  modulus.  These  are  shown  in  Table  2. 
These  values  were  chosen  for  accuracy,  uniformity,  and  from  an  initial  study 
that  c/b  <  1/2  was  the  most  critical  range.  Concerning  accuracy,  it  was  seen 
above  that  c/b  *  1/2  produced  a  maximum  c  for  the  16  inch  base  of  8  inches 
and  a  minimum  e of  0.03125,  which  is  very  near  the  limit  for  accuracy  of 

0.02.  Actually,  the  absolute  greatest  value  of  c  possible  is  16  inches  in 
the  16  inch  base,  and  for  the  element  of  Figure  6  this  corresponds  to  ant 
of  0.015625,  which  is  quite  close  to  0.02.  Therefore,  accuracy  is  not  a 

compelling  reason  to  keep  c/b  less  than  1/2,  although  for  c/b  greater  than 

1/2  one  approaches  the  ragged  edge  of  accuracy. 

The  uniformity  criteria  is  mainly  for  convenience.  For  purposes  of 

comparing  trends  in  Kj  distributions  from  one  base-overlay  combination  to 
another,  it  is  very  helpful  to  have  results  computed  for  the  same  or  similar 
c/b  values.  The  differences  seen  in  Table  2  for  the  12  inch  base  are  because 
of  ease  in  mesh  creation. 

The  main  reason  why  c/b  <  1/2  values  are  used  is  that  is  was  thought  that 
they  were  the  most  critical  values,  especially  those  <  1/6.  This  thought 
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Values  Used  in  Analylsis 

_ C/B  VALUES _ 

1/12,  1/8,  1/6,  7/16,  1/2 
1/12,  5/36,  1/6,  5/12,  1/2 


3] 


CHAPTER  III:  SOLUTION  PROCEDURE 


Development  of  Stress  Distributions 

The  first  step  in  the  solution  process  is  the  calculation  of  the  stress 
distributions  to  be  used  in  superposition  phases  discussed  in  Chapter  I.  In 
this  problem  there  is  no  crack,  and  thus  the  finite  element  program  with  no 
crack  element  can  be  used  to  calculate  the  stress  distribution.  A 
distribution  must  be  found  for  each  combination  of  base  thickness,  overlay 
thickness,  and  modulus  to  be  analyzed.  The  following  combinations  are 
analyzed  in  this  study:  all  modulus  conditions  for  all  base  thicknesses  with 
no  surface,  and  all  modulus  conditions  for  a  base  thickness  of  8  inches  and 
surfacing  of  3  and  6  inches.  This  information  is  summarized  in  Table  3. 

These  specific  combinations  are  analyzed  so  that  all  base  thicknesses  for  a 
given  surface  (0  inches)  can  be  compared,  and  all  surfaces  for  a  given  base 
thickness  (8  inches)  can  be  compared.  A  sample  mesh  is  shown  in  Figure  10 
for  reference.  In  this  mesh,  the  base  thickness  is  equal  to  8  inches  and  the 
overlay  thickness  is  3  inches.  The  input  corresponding  to  this  mesh  is 
listed  in  Appendix  B,  along  with  the  corresponding  output  file.  Note  that 
the  different  modulus  conditions  can  be  achieved  by  simply  changing  the  input 
values.  Also  note  that  the  thickness  of  the  meshes  in  this  section  is  one 
inch,  allowing  the  270  psi  load  to  be  applied  as  a  traction  of  the  same 
magnitude . 

The  results  are  shown  in  Table  U.  The  numbers  in  the  headings  are 
element  numbers,  and  correspond  to  the  elements  in  the  base,  along  the  line 
of  sjmimetry,  where  the  cracks  would  run.  Lower  element  numbers  refer  to 
elements  at  the  bottom  of  the  base,  while  higher  numbers  refer  to  elements  at 
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Table  4.  Stress 

01 strlbutions  for 

Superposition  Phase  one,  psi 

8  Inch 

base,  1 

^0  overlay 

MC 

266 

285 

304 

323 

1 

1705.1 

546.4 

-584.5 

-1772.1 

2 

1229.7 

383.75 

-438.09 

-1314.2 

3 

1455.1 

461.01 

-507.06 

-1530.6 

4 

1032.3 

315.74 

-378.84 

-1126.6 

5 

1086.8 

334.56 

-395.06 

-1178.2 

6 

751.6 

218.24 

-297.81 

-864.87 

12  Inch  base,  no  overlay 


MC 

266 

285 

304 

323 

342 

361 

1 

1205. 5 

704.01 

234.94 

-230.70 

-725.58 

-1279.1 

2 

942.39 

545.64 

178.77 

-185.68 

-578.79 

-1028.9 

3 

1078.3 

627.40 

207.79 

-208.84 

-654.43 

-1157.7 

4 

807.06 

464.26 

149.79 

-162.89 

-503.95 

-900.41 

5 

845.81 

487.56 

158.10 

-169.37 

-525.31 

-936.99 

6 

597.49 

338.45 

104.64 

-128.65 

-389.89 

-704.47 

MC 

266 

285 

16  Inch 
304 

i  base, 

323 

no  overlay 
342 

361 

380 

399 

1 

825.67 

572.11 

340.62 

120.46 

-98.98 

-331.43 

-597.65 

-915.71 

2 

699.45 

482.34 

285.93 

99.95 

-85.69 

-284.17 

-515.73 

-798.05 

3 

786.64 

531.53 

315.90 

111.20 

-92.96 

-310.04 

-560.58 

'862.48 

4 

620.79 

426.45 

251.89 

87.61 

-77.50 

-254.89 

-464.93 

-725.03 

5 

644.37 

443.19 

262.09 

91.00 

-79.94 

-263.64 

-480.13 

'746.89 

6 

478.20 

325.34 

190.36 

63.95 

-62.95 

-202.42 

-373.73 

-593.71 

TABLE  1*2(0):  8  Inch  base,  3  Inch  overlay 


MC 

266 

285 

304 

323 

1 

1466.0 

629.89 

-172.65 

-998.51 

2 

1075.8 

457.60 

-130.19 

-741.01 

3 

1232.9 

666.90 

-66.59 

-717.41 

4 

887.08 

402.36 

-53.30 

-526.53 

5 

903.29 

469.46 

69.95 

-345.38 

6 

632.18 

321.91 

37.94 

-255.06 

8  Inch  base,  6  Inch  overlay 
WC  266 _ 285 _ 304 _ 323 

1  1190.6  615.34  69.38  -382.55 

2  911.27  467.40  50.40  -308.17 

3  1003.4  561.06  146.45  -376.58 

4  751.0  415.71  105.63  -296.51 

5  745.18  468.70  216.88  -312.75 

6  535.85  332.48  150.93  -236.52 


37 

Che  Cop  of  Che  base.  The  numbers  below  Che  heading  MC  refer  Co  Che  modulus 
condicion,  all  of  which  are  depicCed  previously  in  Table  1.  The  numbers 
below  Che  elemenc  headings  are  Che  sCresses  aC  Che  elemenC  cenCroids,  in 
psi.  A  posicive  value  denoces  Cension.  The  sCresses  for  modulus  condicion 
one  in  each  Cable  are  graphed  in  Figure  11. 

Linearize  SCress  DisCribuCions 

As  can  be  seen  from  Figure  11,  Che  scress  disCribucions  are  very  nearly 
linear,  especially  in  Che  range  of  Y  coordinaCe  values  which  encompass  Che 
crack  lengChs  Co  be  invescigaced ,  up  Co  1/2  Che  base  chicknesses.  As  a 
resulc,  for  Che  purpose  of  applying  Che  reverse  sCresses  Co  Che  crack 
lengChs,  ic  is  expedienc  Co  linearize  Che  sCress  disCribucions  by  assuming  a 
linear  discribucion  beCween  each  elemenC  cencroid.  Noce  Chac  elemenc 


cencroids  are  Cwo  inches  aparC  in  Figures  10  and  11.  This  is  because  each 
horizoncal  line  of  elemenCs  in  Che  base  is  Cwo  inches  chick,  for  all  meshes. 
Figure  10  shows  an  8  inch  base.  Addicional  elemenCs  for  12  and  16  inch  bases 
shown  in  Table  A  reflecC  Chese  supplemencal  lines  of  elemenCs. 

The  linearized  disCribucions  are  given  in  Table  5.  LisCed  are  Che 
slopes,  m,  and  Che  inicial  sCress  value  for  each  incerval  beCween  elemenc 
cencroids.  Values  are  lisCed  for  each  modulus  condicion.  The  inCervals  are 
delineaCed  by  che  Y  coordinaCe  values  of  Che  appropriace  elemenc  cencroids. 
Only  inCervals  up  Co  1/2  Che  base  Chickness  are  lisCed.  Noce  Chac  Che 
inCerval  beCween  Che  lowesC  elemenC  cenCroids,  ChaC  is  beCween  51  and  53 
inches,  is  excended  co  Che  boccom  of  Che  base,  which  is  locaced  ac  a  Y  value 
of  50  inches.  The  sCress  aC  Che  bocCom  of  Che  base  is  unknown,  and  Che  same 
slope  chac  aces  beCween  51  and  53  inches  is  assumed  Co  hold  beCween  50  and 
51.  The  scress  ac  che  boCCom  of  Che  base,  or  che  value  for  che  inicial 


Tabic  5.  Linctrixcd  Stress  Distributions,  psi 

8  Inch  base,  no  overlay 
50-53  Inches  53-55  Inches 


579.35 

2284.5 

565.45 

546.40 

422.98 

1672.7 

410.92 

383.75 

497.05 

1952.1 

484.04 

461.01 

358.28 

1390.6 

347.29 

315.74 

376.12 

1462.9 

364.81 

334.56 

266.68 

1018.3 

258.03 

218.24 

12  inch  base,  no  overlay 
50-53  Inches  53-55  Inches  55-57  Inches 


250.75  1456.3 
198.38  1140.8 
225.45  1303.8 
171.40  978.50 
179.13  1024.9 
129.52  727.01 


1456.3  234.54  704.01  232.82  234.94 

1140.8  183.44  545.64  182.23  178.77 

1303.8  209.81  627.40  208.32  207.79 

978.50  157.24  464.26  156.34  149.79 

1024.9  164.23  487.56  163.74  158.10 


338.45  116.65 


16  Inch  base,  no  overlay 

50-53  Inches  53-55  Inches  55-57  Inches  57-59  Inches 


8  Inch  base.  3  Inch  overlay 
50-53  Inches  53-55  Inches 


1 

418.06 

1884.1 

401.27 

629.89 

'  2 

309.10 

1384.9 

293.90 

457.60 

’  3 

333.0 

1565.9 

316.75 

566.90 

1  4 

242.36 

1129.4 

227.83 

402.36 

r  5 

216.92 

1120.2 

201.76 

469.46 

‘ 

155.14 

787.30 

141.99 

321.91 

8  Inch  base.  6  inch  overlay 
50-53  Inches  53-55  Inches 


1 

287.63 

1478.2 

272.98 

615. 

34 

2 

221.94 

1133.2 

208.50 

467. 

40 

3 

221.17 

1224.6 

207.31 

561. 

06 

4 

167.65 

918.65 

155.04 

415. 

71 

1  5 

138.24 

883.42 

125.91 

468. 

70 

1 

a 

101.69 

637.54 

90.78 

332. 
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interval  of  50-53  inches,  is  obtained  by  adding  this  slope  to  the  stress  at 
the  first  element  centroid,  located  a  unit  distance  away  at  51  inches.  Since 
there  is  a  material  (and  thus  stress)  discontinuity  between  the  base  and  the 
subgrade  it  is  impossible  to  obtain  the  stress  at  the  bottom  of  the  base  from 
the  nearest  subgrade  element.  It  is  thus  necessary  to  assume  something,  and 
the  above  procedure  is  consistent  with  the  concept  of  linearized  stress 
distr ibut ions . 

Mesh  Considerations  for  Finite  Element  Runs 

The  sample  mesh  shown  in  Figure  12  is  for  a  base  of  8  inches,  a  surface 
of  3  inches,  and  a  crack  length  of  3.5  inches  (  c/b  of  0.4375).  Also  shown 
is  an  enlargement  of  a  section  of  mesh  in  the  base  which  is  to  fine  to 
represent  with  the  rest  and  an  enlargement  of  the  area  immediately 
surrounding  the  crack  element.  This  crack  element  and  the  area  immediately 
surrounding  it  is  common  to  all  meshes.  There  is  a  mesh  for  each  crack 
length  in  each  base-overlay  combination.  The  common  crack  element  enables 
one  to  build  one  mesh  from  another  quickly  and  easily  and  to  be  able  to 
efficiently  execute  the  analysis.  As  discussed  in  Chapter  II,  there  are  five 
crack  lengths  to  be  investigated  for  each  base-overlay  combination,  and  there 
are  five  base-overlay  combinations  targeted  for  analysis.  Thus,  there  are  25 
separate  meshes  to  be  considered  and  the  efficiency  afforded  by  a  common 
crack  element  is  essential. 

The  maximum  aspect  ratio  in  these  meshes  is  40:1.  This  occurs  as  far 
away  from  the  crack  element  as  possible.  Near  the  crack,  the  aspect  ratios 
are  very  nearly  1:1.  Thus,  near  the  crack  where  greater  accuracy  is 
required,  the  mesh  is  in  its  most  accurate  configuration,  and  far  from  the 
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crack  where  accuracy  is  not  critical,  the  mesh  is  in  its  least  accurate 
configuration. 

The  thickness  of  all  the  meshes  of  this  section  is  one  inch.  Therefore, 
all  stresses  applied  as  loads  transfer  as  tractions  of  the  same  magnitude. 

The  input  file  corresponding  to  Figure  11  is  listed  in  Appendix  B. 

Applying  Reverse  Stresses  to  Crack  Lengths 

The  linearized  stresses  shown  above  are  applied  to  the  meshes  created 
above  with  a  reverse  sign.  Note  that  there  is  one  mesh  for  each  base-surface 
combination  and  one  linearized  stress  distribution  for  each  modulus  condition 
of  each  base-surface  combination.  Therefore,  there  are  25*6  *  150  different 
files  to  be  run  in  the  finite  element  program.  Reverse  stress  values 
calculated  for  the  mesh  shown  in  Figures  11  with  modulus  condition  one  are 
given  in  Table  6.  These  values  also  appear  in  the  input  file  for  the  mesh  in 
Figure  11  that  is  listed  in  Appendix  B.  Note  that  the  values  in  the  input 
file  have  units  of  force/length,  not  stress,  but  have  the  same  magnitude  as 
the  stresses  in  Table  6.  This  is  because  the  thicknesses  of  all  meshes  is 
one  inch.  A  positive  sign  indicates  a  stress  in  the  positive  X  direction. 

Uncorrected  Kj  Values 

Each  of  the  150  runs  took  between  1.50  and  3.00  minutes  on  a  Harris  800 
super  mini-computer.  The  results  are  listed  in  Table  7.  The  values  in  the 
tables  are  uncorrected  stress  intensity  factors  in  ps V*^inch .  The  negative 
signs  indicate  that  the  reverse  stresses  act  to  close  the  crack.  In 
actuality,  these  stress  intensity  factors  are  positive,  because  the  crack  is 
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Table  7.  Uncorrected  Stress  Intensity  Factors.  psI/YTn 
8  inch  base,  no  surface 
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T7T 


MC 

1/12 

1/8 

1/6 

7/16 

1/2 

1 

-916.20 

-1263.99 

-1490.41 

-2372.08 

-2387.92 

2 

-653.28 

-891.83 

-1044.41 

-1585.38 

-1576.26 

3 

-778.18 

-1068.65 

-1256.70 

-1959.78 

-1962.35 

4 

-543.48 

-735.67 

-857.06 

-1257.04 

-1238.55 

5 

-573.83 

-778.89 

-908.91 

-1347.53 

-1331.49 

6 

-387.04 

-512.59 

-589.87 

-801.68 

-773.87 

12  inch 

base,  no  surface 

C/B 

MC 

1/12 

5/36 

1/6 

5/12 

1/2 

1 

-798.13 

-1147.78 

-1216.62 

-1490.31 

-1388.27 

2 

-613.84 

-878.07 

-927.64 

-1105.32 

-1017.78 

3 

-709.19 

-1017.63 

-1077.19 

-1304.38 

-1209.25 

4 

-518.45 

-738.42 

-778.06 

-906.79 

-827.16 

5 

-545.82 

-778.48 

-820.98 

-963.63 

-881.70 

6 

-369.66 

-520.76 

-545.04 

-601.43 

-535.64 

16  inch 

base,  no  surface 

-  T/i - 

4C 

1/12 

1/8 

1/6 

7/16 

1/2 

1 

-603.00 

-779.24 

-913.42 

-870.49 

-783.82 

2 

-502.45 

-647.50 

-757.11 

-704.63 

-631.38 

3 

-572.40 

-736.91 

-860.32 

-800.74 

-725.99 

4 

-439.56 

-565.24 

-659.70 

-602.85 

-537.92 

5 

-458.43 

-589.91 

-688.-91 

-633.19 

-535.64 

6 

-325.38 

-416.05 

-483.27 

-420.69 

-371.11 

inch  base.  3  inch  surface 


MC 

1/12 

1/8 

1/6 

7/16 

1/2 

1 

-763.07 

-1165.89 

-2161.30 

-1962.77 

-2009.58 

2 

-553.75 

-837.89 

-1540.76 

-1360.62 

-1382.10 

3 

-630.24 

-957.27 

-1758.64 

-1593.91 

-1633.61 

4 

-446.75 

-669.68 

-1218.07 

-1066.60 

-1082.68 

5 

-443.83 

-664.67 

-1198.61 

-1083.01 

-1110.89 

6 

-303.38 

-444.61 

-791.12 

-682.53 

-691.41 

8  inch  base.  6  inch  surface 


MC 

1/12 

1/8 

1/6 

7/16 

1/2 

1 

-596.74 

-906.85 

-1655.76 

-1532.27 

-1579.01 

2 

-452.21 

-681.49 

-1237.18 

-1122.00 

-1150.22 

3 

-490.88 

-741.00 

-1338.24 

-1243.83 

-1284.41 

4 

-362.56 

-541.13 

-970.67 

-880.61 

-903.77 

5 

-348.11 

-517.56 

-915.72 

-855.82 

-884.99 

6 

-245.00 

-357.48 

-627.08 

-567.97 

-582.97 
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opened  by  the  loads  on  the  pavement.  However,  as  mentioned  previously,  the 
magnitudes  of  these  numbers  are  correct. 

A  sample  output  file  is  given  in  Appendix  B.  This  output  file 
corresponds  to  the  input  file  in  Appendix  B  that  was  discussed  previously. 

Stress  Intensity  Correction  Factors 

The  equation  used  to  compute  the  correction  factors  is  given  by  equation 

(2)  of  Chapter  I.  In  this  equation,  is  the  stress  in  the  vicinity  of  the 

crack  tip  that  would  be  applied  to  the  portion  of  the  crack  within  the  crack 

element,  if  that  were  allowed.  Because  the  reverse  stress  distribution  is 

assumed  to  be  linear  with  depth,  there  is  no  one  value  of  reverse  stress  in 

the  vicinity  of  the  crack  tip.  is  then  defined  as  the  average  value  of 

the  reverse  stress  that  would  be  applied  to  the  length  of  the  crack  within 

the  crack  element.  The  average  reverse  stress  is  the  stress  that  would  be 

applied  halfway  from  the  crack  tip  to  the  edge  of  the  crack  element.  If  Y  is 

the  vertical  coordinate  of  the  crack  tip,  then  a  acts  at  Y  =  Y  -zl2.  The 

e 

value  of  z,  the  length  of  crack  within  the  crack  element,  is  1/4  inches. 

A  negative  correction  factor  indicates  that  the  reverse  stresses  to  be 
applied  act  in  a  manner  which  closes  the  crack.  They  add  with  the  negative 
stress  intensity  factors  computed  in  the  previous  section  to  increase  the 
magnitude  of  Kj.  Positive  correction  factors  add  to  decrease  fhe  magnitude 
of  the  stress  intensity  factor. 

The  computed  correction  factors  are  given  in  Table  8.  The  computed 
stress  intensity  factors  are  added  to  the  correction  factors  in  order  to  find 
the  final  values  of  Kt. 


The  final  values  of  Kj  are  given  in  Table  9  and  illustrated  in  Figure 
13  through  Figure  17.  Note  that  the  values  of  Kj  are  positive.  The  actual 
stress  intensity  factors  are  equal  in  magnitude  but  opposite  in  sign  to  the 
computed  ones;  that  is,  they  indicate  that  the  crack  will  open  under  load, 
not  close. 

From  the  tables  and  graphs  it  can  be  seen  that  the  Kj  values  for  the  8 
inch  base,  no  surface,  are  greater  than  those  of  the  other  pavements  with  no 
surfaces.  The  16  inch  base  has  the  least  Kj  values  of  the  above 
pavements.  This  is  intuitively  reasonable.  Adding  a  3  inch  surface  to  the  8 
inch  base  generally  decreases  the  stress  intensity  factors,  although  the  Kj 
values  for  c/b  =  1/6  are  greater.  For  the  6  inch  surface,  the  Kj  values 
are  less  than  their  counterparts  in  the  no  surface  case  for  all  the  c/b 
values  for  which  stress  intensity  factors  are  calculated. 


CHAPTER  IV:  CALCULATING  FRACTURE  LIFE 


Regression  Equations  for  Stress  Intensity  Factors 

In  order  to  compute  a  fracture  life  for  the  pavement  system  the  stress 
intensity  factor  distribution  must  be  calculated  for  every  combination  of 
pavement  and  loading,  and  then  this  distribution  must  be  used  in  the  Paris 
equation.  A  simpler  approach  is  to  develop  a  series  of  regression  equations 
to  describe  the  distribution  of  the  stress  intensity  factor.  This  eliminates 
the  need  to  rerun  the  finite  element  program  a  large  number  of  times  to 
develop  a  series  of  stress  intensity  factors  differing  only  slightly  from 
each  other.  The  equation  must  relate  Kj  to  c/b. 

Examining  the  graphs  of  the  stress  intensity  factor  distribution  shown  in 
the  previous  chapter,  it  is  reasonable  to  assume  a  cubic  distribution  of 
versus  c/b,  especially  for  the  cases  with  no  overlay.  In  order  to  define 
this  assumed  cubic  relationship,  regression  techniques  are  used.  A 
statistical  package  on  the  IBM  personal  computer  was  utilized  to  find  the 


best  fit  coefficients 

to 

a  cubic  distribution. 

The 

results  of 

the  polynomial 

regression 

analy 

sis 

are 

given 

in  Table 

10. 

The 

r2  parameters  of 

the 

regression 

are 

also 

listed . 

A  value 

of 

r2  -  1 

.0  denotes 

a  perfect 

fit. 

Note  that 

the 

fit 

is 

very 

good  for 

the 

cases 

with  no 

surface , 

and 

satisfactory  for  those  with  a  surface.  Plots  of  Kj  versus  c/b  for  both  the 
cubic  distribution  and  the  distributions  of  the  previous  section  are  shown 
for  representative  cases  in  Figure  18  through  Figure  22. 

For  reasons  discussed  previously,  the  range  of  c/b  values  for  which 
stress  intensity  factors  were  computed  extend  from  1/12  to  1/2.  For  no  c/b 
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v  v.  -r.  ' 


Table  10.  Regression  Coefficients  for  Cubic  Polynomial 
Fit  to  Stress  Intensity  Distribution. 

K  •  B  ♦  B  (C/B)  ♦  B  (C/B)  ♦  B  (C/B) 


2047.043 

1490.105 

1753.576 

1260.895 

1351.920 

938.9210 


8  Inch  base,  no  surface 

-ilZI  B?  Bt 


6567.108 

4544.518 

5513.253 

3666.674 

3600.485 

2345.091 


-14172.16 

-10815.35 

-12415.28 

-9340.74 

-8792.111 

-6942.935 


5067.807 

4237.885 

4625.061 

3857.532 

3060.146 

3093.236 


0.9861128 

0.9940159 

0.9903571 

0.9966238 

0.9971226 

0.9988023 


1297.536 

1019.034 

1162.857 

875.8471 

916.8071 

654.5362 


12  Inch  base,  no  surface 
ilH  B2_  B;} 


7693.483 

5847.388 

6804.254 

4887.160 

5163.552 

3373.223 


-22383.32 

-17613.80 

-20087.00 

-15129.88 

-15845.97 

-11169.04 


14619.22 

11676.70 

13202.72 

10145.02 

10587.44 

7681.451 


0.9951294 

0.9962779 

0.9957248 

0.9968587 

0.9966850 

0.9978152 


16  inch  base,  no  surface 
Bi  B?  Bt _ 


808.9672 

687.1326 

776.4540 

611.02611 

633.80810 

472.87540 


6928.043 

5742.700 

6547.792 

5000.017 

5223.816 

3650.309 


-22642.82 
^19189.89 
•22345.70 
-16980. B5 
-17649.63 
-12930.04 


17098.41 

14661.88 

17517.90 

13073.25 

13556.44 

10141.79 


0.9999410 

0.9998762 

0.9998906 

0.9999341 

0.9999398 

0.9998798 


•394.5381 

-235.5056 

-280.6712 

-137.2445 

-112.6183 

-12.82105 


8  Inch  base. 

jir — z! 


3  Inch  surface 

[iz  zim 


37465.05 

26687.00 

30515.79 

20948.11 

20547.46 

13339.98 


-117437.1 

•84582.00 

•96126.81 

-66598.86 

-64805.32 

-42658.10 


105884.90 

76535.600 

87240.120 

60381.990 

59002.900 

38802.090 


0.8500299 

0.8609347 

0.8503666 

0.8641640 

0.8512328 

0.8696607 


-247.4966 

-148.9915 

-153.3865 

-72.99951 

-35.03687 

17.19202 


8  Inch  base. 

jiz — z! 


6  Inch  surface 

EiZ  bt 


28694.99 

21328.18 

23050.09 

16591.39 

15331.10 

10486.71 


-89707.74 

-67057.85 

-71954.05 

-52166.28 

-48342.45 

-32985.49 


81544.99 

60936.81 

65549.56 

47507.98 

44153.58 

30127.45 


0.8456994 

0.8496895 

0.8465933 

0.8494838 

0.8497349 

0.8487314 


V  'V 
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Figure  18.  Comparison  of  Computed  and  Regression  Stress  Intensity  Factors 


Distribution  of  Stress  Intensity  Factor 
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values  greater  than  1/2  were  stress  intensity  factors  calculated.  Thus,  for 
c/b  greater  than  1/2,  the  assumed  cubic  distributions  give  the  only  Kj 
values;  they  must  be  relied  on  exclusively.  The  graphs  of  the  predicted 
distributions  for  c/b  greater  than  1/2  are  shown  in  Figure  23  through  27.  As 
can  be  seen,  the  plots  for  all  cases  get  smaller  directly  after  c/b  =  1/2, 
and  then  double  back  and  increase.  For  the  8  inch  base,  no  surface,  the  Kj 
values  go  negative  and  never  become  positive  again.  For  the  rest  of  the 
cases,  Kj  is  always  positive  past  c/b  =  1/2,  and  for  the  8  inch  base,  3  and 
6  inch  surface  cases  Kj  gets  very  large  at  c/b  =1.  In  addition,  noting 
that  Bq  equals  Kj  at  c/b  =0,  it  can  be  seen  that  for  the  8  inch  base,  3 
and  6  inch  surface  cases  Kj  is  negative  at  c/b  =  0.  The  exception  is 
modulus  condition  6  of  the  6  inch  surface,  where  Kj  manages  to  remain 
positive  at  c/b  =  0. 

The  predicted  Kj  values  at  c/b  =  0  cannot  be  checked  (there  is  no  crack 
length  on  which  to  apply  reverse  stresses),  but  for  c/b  greater  than  1/2  the 
cubic  values  should  be  checked  with  calculated  Kj  values.  It  may  be  that 
the  cubic  distributions  are  misleading  in  this  range,  but  these  values  have 
not  been  checked. 

Computation  of  Crack  Propagation  Histories 

To  calculate  the  rate  of  crack  propagation,  the  Paris  equation  described 
in  Chapter  II  was  coded  into  a  Fortran  program  to  allow  easy  application  of 
stress  intensity  factors,  pavement  geometry  and  structural  properties  into  a 
program  that  predicted  the  number  of  load  cycles  to  failure.  Up  to  this 
point,  all  the  necessary  input  has  been  calculated  except  for  the  initial  and 
final  crack  lengths  for  each  modulus  condition,  base,  and  surface 


Figure  Ik.  Distribution  of  Stress  Intensity  Factors  Computed  by 
Regression  Equation  for  Various  Modulus  Conditions 


16  Inch  Base,  No  Surface 


Figure  27.  Distribution  of  Stress  Intensity  Factors  Computed  by  Regression 
Equation  for  Various  Moduli  Conditions 
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combination.  These  are  found  by  calculating  the  roots  of  the  Kj  cubic 
distributions.  The  first  two  roots  in  the  range  0  ^  c/b  £  1  which  define  a 
positive  region  of  are  the  endpoints,  as  explained  in  Chapter  1.  If 
Kx  is  positive  in  the  range  0  ^  c/b  <  1,  then  the  endpoints  are  0  and  1. 
The  endpoints  found  by  using  this  procedure  are  shown  in  Table  11.  When 
input,  if  the  endpoint  in  question  is  a  root,  then  it  is  modified  slightly  so 
that  it  does  not  assign  a  zero  (or  negative  due  to  round-off  errors)  value  to 
Kx  and  thus  cause  computational  errors.  If  it  is  a  right  endpoint,  it  is 
made  smaller,  and  if  it  is  a  left  endpoint  it  is  made  larger  so  that  Kx  is 
always  positive. 

The  crack  propagation  histories  for  modulus  condition  4,  all  base  and 
surface  combinations,  are  given  in  Figure  28  through  Figure  32.  These 
histories  correspond  to  arbitrarily  selected  values  of  A  and  n  equal  to 
10"^®  and  2.0  respectively.  The  numbers  defining  these  graphs  can  be 
found  in  the  output  in  Appendix  B,  as  mentioned  before.  The  histories 
graphed  are  typical  in  their  relationship  to  each  other  for  each  modulus 
condition  and  A  and  n  combination. 

It  appears  from  the  graphs  that  of  the  pavements  with  no  surface,  the  8 
inch  base  is  the  most  critical,  followed  by  the  12  and  then  16  inch  bases. 
Although  the  number  of  cycles  required  to  traverse  the  last  crack  increment 
in  the  8  inch  base  is  five  orders  of  magnitude  greater  than  the  total  to 
failure  of  the  12  and  16  inch  bases,  the  load  cycles  needed  to  reach  near 
failure  in  the  8  inch  base  is  significantly  smaller  than  the  aforementioned 
totals.  It  is  reasonable  to  believe  that  in  such  a  state  of  near  failure 
some  perturbation  could  advance  the  crack  to  failure.  Besides,  the  8  inch 
base  is  practically  destroyed  at  near  failure  anyway. 


t   •  —  ^  • 


Figure  28.  Progression  of  Crack  (c)  Through  Base  (b) 
as  a  Function  of  Load  Cycles,  N 


Cycles,  N  xlO 


c/b 


Figure  31.  Progression  of  Crack  (c)  Through  Base  (b) 
as  a  Function  of  Load  Cycles,  N 


^ 


Figure  32.  Progression  of  Crack  (c)  Through  Base  (b) 
as  a  Function  of  Load  Cycles,  N 


The  3  and  6  inch  surfaces  are  added  to  the  8  inch  base  with  the  intent  of 
increasing  the  fatigue  life.  From  the  graphs  it  is  seen  that  theoretically 
they  do,  but  perhaps  realistically  they  won't.  In  both  the  3  and  6  inch 
surface  cases,  it  takes  a  tremendous  amount  of  load  cycles  to  advance  the 
crack  through  the  first  increment  of  approximately  0.05  inches.  After  that, 
though,  it  takes  less  cycles  to  propagate  the  crack  to  failure  than  it  does 
for  the  crack  to  reach  near  failure  in  the  no  surface  case.  Therefore,  if 
the  crack  were  somehow  jumped  past  the  initial  increment  in  the  3  and  6  inch 
surface  cases,  due  to  say  a  starter  flaw,  which  may  typically  have  dimensions 
capable  of  doing  this,  those  pavements  would  be  more  critical  than  the  no 
surface  case.  It  is  certainly  reasonable  to  think  that  the  first  increment 
could  be  jumped,  as  it  is  assumed  that  a  zone  of  negative  is  jumped  to 
get  to  the  first  increment  in  the  first  place.  In  this  case  the  negative 
zone  is  only  about  0.01  inches  at  most,  so  that  it  is  no  great  achievement  to 
bypass  it.  So,  in  one  sense,  adding  a  surface  will  greatly  postpone  failure 
as  long  as  the  negative  Kj  zone  and  the  first  crack  increment  are  not 
vaulted.  In  another  sense,  if  the  above  were  circumvented,  adding  a  surface 
could  hasten  failure. 

As  would  be  expected,  the  3  inch  surface  case  requires  more  cycles  to 
failure  than  the  6  inch  surface  case  for  the  same  pavement  structure. 

Calculations  With  Actual  Laboratory  Fracture  Properties 

Extensive  laboratory  studies  conducted  at  Texas  A&M  University  have 
indicated  that  the  fracture  parameters,  A  and  n  can  be  predicted  from  tests 
along  with  Kj^.  for  stabilized  materials.  These  parameters  are  necessary  to 
predict  the  fracture  life  of  a  material  insitu.  Typical  values  determined 


are  as  follows: 
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CEMENT  LEVEL 


E=355,600  psi  v  =  0,15 


CREEP  VALUES 


A=8, 12x10' 


Ktp™1 38 . 6 


FATIGUE  VALUES 


n-25.28 


A-15.4x10"27 


n*10.0 


E=597,600  psi  v  *  0.15 


CREEP  VALUES 


A=l. 45x10' 


KTr=209.3 


n=«21 .24 


FATIGUE  VALUES 


A»10xl0' 


n=22.0 


Comparison  of  these  critical  stress  intensity  factors  indicates  that  the 
stress  intensity  factors  calculated  using  the  procedures  outlined  here  are 
all  larger  than  the  critical  stress  intensity  factors  which  means  that 
fracture  will  proceed  instantaneously.  This  indicates  the  inaccuracy  of  the 
pavement  structure  selected  for  modelling  the  stress  intensity  factors 
calculated . 


CHAPTER  V:  CONCLUSIONS  AND  RECOMMENDATIONS 


The  analysis  procedure  presented  in  this  report  illustrates  the 
complexities  involved  in  calculating  stress  intensity  factors,  and  applying 
them  to  a  structure  to  predict  cycles  to  failure.  The  stress  intensity 
factors  used  in  this  analysis  are  not  as  accurate  as  may  be  desired,  due  to 
the  use  of  plane  strain,  plane  stress  solutions  used  in  the  finite  element 
procedure  to  calculate  the  initial  stresses.  The  modification  of  more 
sophisticated  programs  to  accept  the  fracture  elements  in  an  axi-symmetr ic 
analysis  would  increase  the  accuracy  of  the  predictions  of  fracture  in  the 
cement  stabilized  base  materials,  or  any  material  susceptible  to  fracture. 

The  methodology  of  applying  the  principles  of  fracture,  as  demonstrated 
in  the  development  of  this  report  are  applicable  to  a  design  methodology,  and 
indicate  the  modularized  approach  which  must  be  taken  to  adequately 
characterize  the  pavements,  their  geometry,  and  the  loading  conditions  which 
may  be  expected  to  occur  in  an  actual  airfield  pavement. 

The  principles  of  fracture  mechanics  and  crack  propagation  can  be  applied 
to  pavement  design  life,  to  account  for  the  gradual  failure  due  to  repeated 
loading  over  fracture  susceptible  materials.  The  parameters  required  to 
predict  crack  propagation  are  the  material  properties  of  A  and  n,  as 
described  in  this  report.  A  significant  amount  of  testing  has  gone  into 
evaluating  these  parameters  for  cement  stabilized  materials,  and  the 
influence  of  these  material  properties  on  fracture  life  was  shown. 

The  predicted  stress  intensity  factors  were  high,  due  directly  to  the 
level  of  structure  modeled  and  the  magnitude  of  loading  investigated. 


A  broader  range  of  loading  conditions,  and  a  thicker  pavement  structure  would 
have  produced  more  moderate  stress  intensity  factors  which  would  have  shown  a 
more  pronounced  effect  on  fracture  life  than  was  seen  in  this  report. 

1.  The  methodology  developed  in  this  study  indicates  the  potential  for 
considering  fracture  in  pavement  design  methodology. 

2.  The  stress  calculation  scheme  can  be  refined  and  replaced  by  a 
general  program  to  analyze  an  axisymmetric  problem  with  multiple 
wheel  loadings.  This  can  be  a  simple  elastic  layer  program,  a 
stress  dependent  elastic  layer  program,  or  even  a  sophisticated 
stress  dependent  finite  element  program  such  as  Illi-  Pave. 

3.  The  hybrid  crack  tip  element  provides  a  simple  means  of  accurately 
evaluating  the  influence  of  loading  on  fracture  life  in  fracture 
susceptible  layers  of  a  pavement  system  given  the  stresses  present 
in  the  layers. 

4.  The  modular  approach  developed  to  provide  stress  intensity  factors 
using  regression  techniques  to  replace  repetitive  calculations  with 
the  finite  element  program  is  an  accurate  procedure  to  obtain  the 
stress  intensity  factors  in  a  manner  applicable  to  design  without 
extensive  use  of  computer  execution  time. 
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APPENDIX  A:  NOTATION  AND  CONVERSION  FACTORS 


MEANINGS 

regression  constant  in  fracture  equation;  area  in  other  cases 
Element  areas  in  finite  element  development 
base  course  thickness  in  fracture  analysis 
crack  length  within  a  pavement  layer 
stress  intensity  correction  factor 
Elastic  Compliance  Tensor 
2.71828 

Young's  modulus  of  elasticity 
shear  modulus 
stress  intensity  factor 
creep  exponent 

regression  constant  (exponent)  in  fracture  equation 

number  of  cycles  in  fracture 

load 

nodal  displacements 

statistical  correlation  coefficient 

radius 

Surface  for  integration 

Traction  vector  normal  to  integration  path 

boundary  displacement 

xj,  Yi  directions  respectively 

error  in  stress  intensity  predictions,  strain 

local  coordinate  axis 

angle 

Poisson's  ratio 
stress 

principal  stress 
shear  stress 

summation 

square  root 

natural  (Naperian)  logarithm  (base  e) 

Strain  energy  release  rate 
coordinate  distance  from  crack  tip 
stress  applied  to  crack  tip 


APPENDIX  B :  DESCRIPTION  OF  PROGRAM  UNITS 


Finite  Element  Fracture  Program 


A  flow  diagram  of  subroutines  is  shown  in  Figure  B-1 .  The  following  is  a 
brief  description  of  each  subroutine  and  the  input  data  format  follows. 
DATAIN: 

Reads  and  echo  prints  all  data  pertaining  to  the  conventional 
elements.  It  also  performs  checks  on  this  data. 

ASEMBL; 

Initializes  and  assembles  the  glabal  stiffness  matrix  (including  the 
crack  element)  and  the  global  load  vector.  Modifies  the  above  to 
reflect  the  geometric  boundary  conditions. 

QUAD: 

Computes  the  stress-strain  matrix,  stiffness  matrix,  body  force 
vector,  and  strain-displacement  matrix  of  either  a  4-CST 
quadrilateral  or  a  triangular  element. 

CST: 

Computes  the  strain-displacement  matrix,  stiffness  matrix,  and  body 
force  vector  of  a  CST  element. 


GEOMBC : 


Applies  prescribed  displacement  boundary  conditions  at 
specified  node. 


a  single. 


Reads  the  input  data  pertaining  to  the  crack  element.  Computes  the 


constants  and  of  Section  II.  Incorporates  the  crack  element 
stiffness  matrix  into  the  global  stiffness  matrix  which  had 
previously  contained  stiffnesses  from  only  conventional  elements. 

HYBRID: 

Computes  the  matrices  [Gl  and  [H]  of  Section  II.  Computes  the  crack 
element  stiffness  matrix.  Computes  the  (Bj)^  and  (Bjj)^  row 
vectors  of  Section  II. 

LOC: 

The  crack  element  stiffness  matrix  is  not  derived  in  two-dimensional 
array  form  but  rather  in  vector  form.  Thus,  in  CRACK  and  SIFAC  where 
this  stiffness  matrix  is  used  in  calculations  and  manipulations,  some 
operator  must  be  used  to  translate  from  two-  to  one-dimensional 
form.  This  operator  is  LOC.  'Oiis  is  probably  done  to  save  memory. 

MFSD: 

Uses  numerical  methods  to  perform  the  integration  implied  by 
equations  (16)  and  (17)  of  Section  II. 

SINV: 

Performs  the  inversion  of  [H] . 

BANSOL : 


Triangular izes  the  global  banded  stiffness  matrix  by  symmetric 
Gauss-Doolittle  decomposition  and/or  solves  for  the  global 
displacement  vector  corresponding  to  a  given  load  vector 


STRESS: 


Computes  strains,  stresses,  and  principal  stresses  for  conventional 
elements.  Computes  the  strain  energy  due  to  all  the  conventional 
elements  which  can  be  used  to  calculate  stress  intensity  factors  by 
energy  methods,  as  mentioned  in  Section  I.  Prints  stresses  and 
principal  stresses  at  element  centroids,  and  the  strain  energy. 

SIFAC: 

Computes  the  stress  intensity  factors  by  equations  (23)  in  Section 
II.  Computes  the  crack  element  strain  energy  and  the  global  strain 
energy,  which  can  be  used  to  compute  stress  intensity  factors  by 
energy  methods,  as  in  Section  I.  Prints  the  crack  element  strain 
energy. 

Input  Guide-  Descriptions  are  included  for  each  card. 

IDENTIFICATION  CARD: 

One  card  per  problem:  Format  (l5,3x,9A8) 

cc  1-5  NPROB:  The  problem  number.  If  NPROB  -  0,  execution  of 

the  program  is  halted. 

cc  9-80  TITLE(I):  Title  of  the  problem.  The  vector  has  9  elements, 

each  of  which  contains  9  characters  of  the  title. 

Multiple  problems  can  be  handled,  but  as  soon  as  the  next  problem  is  read 

the  previous  problem  is  lost,  so  a  problem  cannot  be  recalled. 


BASIC  PARAMETERS: 


One  card  per  problem:  Format  615 

cc  1-5  NNP:  Number  of  nodal  points 

cc  6-10  NEL:  Number  of  conventional  elements 

cc  11-15  NMAT:  Number  of  different  materials 

cc  16-20  NLSC:  Number  of  surface  tractions 

cc  21-25  NOPT:  Option  for  stress  state.  1  =  plane  strain,  2  =  plane 

stress . 

cc  26-30  NBODY:  Option  for  body  force.  0  =  no  weight,  1  =  weight  in 

the  negative  y  direction, 
cc  31-35  NCKEL:  Number  of  crack  elements. 

MATERIAL  PROPERTIES 

NMAT  cards  per  problem:  Format  4F10.0 


cc 

1-10 

E: 

Modulus  of  elasticity 

cc 

11-20 

PR: 

Poisson's  ratio 

cc 

21-30 

RO: 

Density  of  the  material 

cc 

31-40 

TH: 

Thickness  of  the  material 

NODAL  POINT  DATA: 

Up  to  NNP  cards  per  problem:  Format  (215,4F10.0) 
cc  1-5  M:  Nodal  point  number 

cc  6-10  KODE(I):  Index  of  displacement  and  concentrated  load 

conditions  at  note  I.  The  values  Kode  can  assume  and 
the  conditions  assigned  to  each  are  given  in  Table 
IV- 1. 

cc  11-20  X:  Horizontal  coordinate  of  node  I. 

cc  31-30  Y:  Vertical  coordinate  of  node  I. 
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cc  31-40  ULX:  Concentrated  load  or  displacement  in  X  direction  at  node 
I. 

cc  41-50  ULY:  Concentrated  load  or  displacement  in  Y  direction  at  node 
I. 

Usually  one  card  is  needed  for  each  node.  However,  if  some  nodes  fall  on 

a  straight  line  and  are  equidistant,  data  for  only  the  first  and  last  points 

of  this  group  are  needed.  Intermediate  nodal  point  data  are  automatically 
generated  by  linear  interpolation.  The  nodal  data  must  be  entered  in  order 
from  smallest  to  largest,  leaving  out  those  nodes  which  are  to  be 
interpolated.  The  nodes  which  are  interpolated  are  assigned  values  of  KODE  = 
0,  ULX  =  0,  and  ULY  *  0.  The  signs  of  presecribed  nodal  displacements  or 
forces  follow  the  signs  of  the  coordinate  directions  assigned  by  the  user. 
ELEMENT  DATA: 

Up  to  NEL  cards  per  problem:  Format  615. 

cc  1-5  M:  The  element  number 

cc  6-10  IE(M,1):  The  index  of  the  first  node  of  a  CST  or  4-CST 

quadrilateral  element. 

cc  11-15  IE(M,2):  The  index  of  the  second  node  of  a  CST  or  4-CST 

quadrilateral  element. 

cc  16-20  IE(M,3):  The  index  of  the  third  node  in  a  CST  or  4-CST 

quadrilateral  element. 

cc  21-25  1E(M,4):  The  index  of  the  fourth  node  in  a  CST  or  4-CST 

quadrilateral  element.  If  it  is  a  CST  element,  the 
index  of  fourth  node  equals  that  of  the  third  node; 
that  is,  IE(M,3)  »  IE(M,4). 

cc  26-30  IE(M,5):  Material  type  number  corresponding  to  element  M. 
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Usually,  one  card  is  needed  for  each  element.  However,  if  some  elements 
are  on  a  line  in  such  a  way  that  their  corner  node  indexes  each  increase  by 
one  compared  to  the  previous  element,  only  the  data  for  the  first  element  on 
the  line  need  be  input.  As  the  elements  on  the  line  are  generated  by  adding 
one  to  each  node  of  the  preceeding  element,  starting  with  the  first  element 
on  the  line,  the  last  element  on  a  line  needn't  be  input.  However,  data  for 
the  last  element  in  the  assemblage  must  be  input  whether  it  could  be 
generated  or  not.  Also,  please  note  that  triangular  elements  cannot  be 
generated  from  quadrilateral  elements  because  the  third  and  fourth  indexes  of 
a  triangular  element  are  equal.  The  same  material  type  as  the  first  element 
on  a  line  is  assigned  to  all  elements  generated  on  that  line. 

For  a  right-handed  coordinate  system,  the  nodal  indices  for  an  element 
must  be  input  counter-clockwise  around  the  element. 

SURFACE  TRACTIONS: 

As  many  cards  as  needed:  Format  (215,4F10.0) 

cc  1-5  ISC:  I,  the  first  node  upon  which  tractions  act. 

cc  6-10  JSC;  J,  the  second  node  upon  which  tractions  act. 

cc  11-20  SURXl;  The  intensity  of  the  traction  in  the  X  direction 

acting  at  node  I. 

cc  21-30  SURX2:  The  intensity  of  the  traction  in  the  X  direction 

acting  at  node  J. 

cc  31-40  SURYl :  The  intensity  of  the  traction  in  the  Y  direction 

acting  at  node  I. 

The  intensity  of  the  traction  in  the  Y  direction 
acting  at  node  J. 


cc  41-50  SURY2: 


Surface  tractions  must  be  specified  between  two  adjacent  nodes  only.  The 
tractions  in  both  the  X  and  Y  directions  are  assumed  to  vary  linearly  between 
the  two  nodes.  Intensities  are  expressed  in  units  of  force/length  so  that 
pressures  must  be  multiplied  by  the  thickness  before  being  input  into  the 
computer.  The  signs  of  the  tractions  follow  the  directions  of  the  coordinate 
axes  assigned  by  the  user. 

CRACK  ELEMENT  DATA: 

There  are  NCKEL  cards  per  problem  for  each  of  the  two  types  of  card  which 
comprise  the  crack  element  data. 

Card  one:  Format  ( 215 , 2F10 .0 , 15 ) 


CC 

1-5 

KEY: 

The  type  of  crack  element.  1  *  five  node  case,  2  = 

nine  node  case. 

cc 

6-10 

MATYP: 

The  type  of  material  where  the  crack  element  is 

placed . 

cc 

11-20 

XC: 

The  horizontal  coordinate  of  the  crack  tip. 

cc 

21-30 

YC: 

The  vertical  coordinate  of  the  crack  tip. 

cc 

31-35 

NCOT: 

Flag  which  determines  which  direction  nodal  indexes 

are  to  be  counted  around  the  crack  element.  For  the 
five  node  case  I  =  clockwise,  0  =  counter-clockwise. 


Card  two:  Format  1015 

cc  (5l-4)-5l  KCRK(l):  Nodal  incidence  for  the  Ith  node  of  the  crack 

element . 

cc  (5K+1 )-5 (K+1 )  MAXDIF:  Maximum  difference  between  the  nodal  incidences 

of  rhe  crack  element.  Used  to  compute  the 
bandwidth  after  the  addition  of  the  crack 


element . 


For  the  five  node  case  there  are  five  incidences  input  for  the  crack 
element  (K  =  5)  and  for  the  nine  node  case  there  are  nine  (K  =  9).  In  either 
case  the  input  for  MAXDIF  immediately  follows  the  last  incidence  input.  All 
crack  element  data  is  omitted  if  NCKEL  =  0. 


The  listing  of  the  finite  element  computer  code  is  given  on  the  following 


‘  fc"** 
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Crack  Propagation  Calculation  Program 


Figure  B-2  shows  the  program  flow  diagram.  The  following  is  a  brief 
description  of  the  main  program  and  each  subroutine. 

CRACK: 

Main  program.  Sets  values  of  necessary  parameters  and  reads  in 
necessary  data.  Iterates  from  problem  to  problem,  reading  in  data 
specific  to  the  problem  and  calling  NVSC.  Writes  out  modulus 
values,  A,  and  n  for  each  modulus  condition  considered. 

NVSC:  Calculates  and  applies  crack  increments.  Computes  Nf  values  for 

each  increment.  Computes  Nf  for  each  crack  increment.  Prints  the 

base  and  overlay  thicknesses  and  the  modulus  condition  of  the 
pavement  being  analyzed. 

TRAPRLE : 

Given  the  left  and  right  limits  of  the  crack  increment,  numerically 
integrates  the  Paris  equation  using  the  trapezoidal  rule,  producing 
Nf.  Iterates  until  the  percent  difference  is  below  10“^  or  to 
a  limit  of  13  iterations. 

SIMPRLE: 

Performs  exactly  as  does  TRAPRLE  except  that  Simpson's  rule  is  used 
for  the  numerical  integration. 


Determines  the  value  of  Kj  for  a  given  value  of  (c/b)  according  to 
equation  (6),  Section  I.  Computes  the  integrand  of  equation  (4), 
Section  I,  including  the  base  thickness  term.  Paris  is  a  function 
which  returns  as  its  value  the  above  integrand. 

Input  Guide-  The  input  parameters  and  formats  are  explained  here. 


BASIC  PARAMETERS: 

One  Card:  Format  215 

cc  1-5  Nmodcon:  The  number  of  modulus  conditions  considered  for  the 

problem. 

cc  6-10  K:  A  value  of  K  ■  0  indicates  that  the  trapezoidal  rule  will 

be  used.  Any  other  5  digit  integer,  normally  1,  indicates 
that  Simpson's  rule  will  be  used. 

Nmodcon  cards:  Format  4F10.0 

cc  1-10  Evalues(l , 1 ) :  Modulus  of  elasticity  of  the  base  corresponding 

to  modulus  condition  1. 

ccll-20  EvaluesC I ,2) :  Modulus  of  elasticity  of  the  overlay 

corresponding  to  modulus  condition  1 
cc  21-30  AA(I):  Values  of  the  parameter  A  in  the  Paris  equation 

corresponding  to  modulus  condition  1. 

Value  of  the  parameter  n  in  the  Paris  equation 
corresponding  to  modulus  condition  1. 


cc  31-40  NN(I): 


PROBLEM  PARAMETERS: 


One  card  per  problem:  Format  315 

cc  1-5  Nprob:  The  problem  number  of  the  current  problem.  If  the 

value  is  0,  execution  of  the  program  is  halted, 
cc  6-10  Base:  Thickness  of  the  base, 

cc  11-15  Olay:  Thickness  of  the  overlay. 

Nmodcon  cards  per  problem:  Format  6F10.0 

cc  1-10  Bq:  Coefficient  Bq  in  equation  (6),  Section  I. 

cc  11-20  Bj :  Coefficient  Bj  in  equation  (6),  Section  I. 

cc  21-30  82:  Coefficient  82  in  equation  (6),  Section  I. 

cc  31-AO  83:  Coefficient  83  in  equation  (6),  Section  I. 

cc  41-50  Left:  Initial  crack  length  for  a  pavement  with  given  base 

and  overlay  thicknesses  and  modulus  condition, 
cc  51-60  Right:  Final  crack  length  for  a  pavement  with  given  base  and 

overlay  thicknesses  and  modulus  condition. 

The  computer  code  listing  for  the  paris  equation  calculation  scheme  is 
given  on  the  following  pages. 
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